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ABSTRACT 


The  research  effort  focuses  on  the  exploration  of  sensing  and  control  techniques 
in  the  framework  of  fusion  and  distributed  control  pertinent  to  the  problem  of  deploying 
autonomous  unmanned  surface  platforms  and  clusters  of  platforms  in  mission 
applications  in  support  of  future  naval  capability.  The  considerations  which  comprise  the 
approach  to  the  distributed  autonomous  problem  are  threefold:  (1)  platform  sensing 
capability  for  environment  and  mission;  (2)  distributed  based  control  of  autonomous 
platforms  and  clusters;  (3)  the  intelligent  fusion  of  sensor  data,  mission  and  control 
action.  Fundamentally,  the  effort  deals  with  the  study  of  pertinent  modeling  approaches 
which  may  be  useful  in  enhancing  the  processing  and  fusion  of  data  to  achieve 
appropriate  distributed  control  action.  This  provides  a  vehicle  for  the  exploration  of 
mission  specific  distributed  autonomous  platform  performance.  From  a  real-time 
application  perspective,  the  real-time  turnaround,  time  correlation  and  the  adaptability  of 
a  fusion  based  distributed  control  are  of  concern.  This  requires  the  development  of  novel 
approaches  to  sensor  return  integration  and  interpretation  in  conjunction  with  novel 
approaches  to  developing  a  fusion  based  distributed  control  strategy  capable  of 
adaptation,  learning  and  information  sharing.  Potential  mission  applications  include  the 
implementation  of  these  techniques  in  unmanned  vehicle  control,  mine  warfare,  anti¬ 
submarine  warfare,  and  anti-terrorism/force  protection. 
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1.0  INTRODUCTION 


The  envisioned  US  Navy  future  naval  capability  incorporates  the  development  of 
autonomous  unmanned  platforms  and  cluster  of  platforms  [5].  This  autonomous  platform 
architecture  may  comprise  of  thousands  of  fully  netted  unmanned  vehicles  and 
unattended  sensors  (e.g.  FORCEnet)  as  illustrated  in  Figure  1 .  Key  to  the  success  of  this 
concept  is  the  ability  to  maintain  a  functional  network  of  autonomous  distributed 
platforms  deployed  as  sensors  and  combatants. 

1.1  Application  Overview 

Crucial  to  the  successful  realization  of  such  a  system  is  the  development  of  the  capability 
of  maintaining  real-time  sensor  data  turnaround  in  concert  with  autonomous  formation 
control  driven  by  mission  requirements  and  environmental  assessment  based  on  sensor 
fusion.  Techniques  applied  to  attain  distributed  autonomous  platform  cluster  survivability 
and  autonomous  formation  control  are  critical  in  accomplishing  a  hostile  environment 
mission. 

The  effort  presented  addresses  these  critical  issues  in  the  context  of  seeking  to 
develop  innovative  approaches  for  sensor  data  fusion  and  in  the  incorporation  of  the 
fused  data  for  adaptive  distributed  autonomous  control  in  formation  control,  addressing 
mission  and  operating  environment.  Pertinent  naval  applications  include  unmanned 
vehicle  control  (e.g.  UV),  mine  and  undersea  warfare  (MUW),  anti-submarine  warfare 
(ASW),  and  anti-terrorism/force  protection  (ATFP). 
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Figure  1:  Illustration  of  large-scale  autonomous  operations. 

1.2  Problem  Description 

The  fundamental  problem  of  investigating  an  approach  to  autonomous  platform 
distributed  formation  control  based  on  aspects  of  sensor  fusion  pertinent  to  platform-to- 
platform  sensing  is  addressed.  The  research  focuses  on  time  delay  aspects  of  sensing  and 
on  real-time  autonomous  intelligent  navigational  trajectory  generation  [19],  [20],  The 
former  is  an  issue  when  the  separation  between  platforms  is  very  large  compared  to  the 
duration  of  signals  generated  in  the  operating  environment.  The  latter  deals  with 
formation  control  of  network  of  autonomous  vehicles  which  are  capable  of  adaptation 
based  on  the  sensed  operating  environment,  condition  of  the  platform  and  mission 
objectives. 

Section  2.0  of  the  report  presents  results  on  autonomous  navigation  and  control. 
The  formation  control  problem  is  investigated  from  the  perspective  of  employing  a  vector 
integration  to  endpoint  differential  structure  to  adaptively  generate  trajectories  while 
maintaining  a  coordinate  map  base  of  the  fusion  environment.  The  combination  of 
sensing,  formation  control  and  mission  fusion  addresses  current  critical  limitations  in 
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implementing  a  full  scale  operational  autonomous  platform  network.  The  modeling  and 
simulation  of  these  processes  will  provide  a  basis  for  performance  assessment  and  future 
development. 

Section  3.0  of  the  report  presents  results  on  the  autonomous  sensing  problem  for 
distributed  agents  with  application  to  autonomous  surface  vehicle  distributed  fluid  supply 
systems.  The  sensing  problem  concerns  itself  with  the  efficient  estimation  of  the  time-of- 
flight  associated  with  acoustic  signals  between  separated  platform  sensors.  The  use  of 
hyper-geometric  properties  will  be  investigated  as  a  means  of  reducing  the  number  of 
sensors  for  identification  of  time-of-flight  and  as  a  means  of  reducing  computational 
complexity  in  fusing  the  sensor-to-sensor  data. 
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2.0  AUTONOMOUS  NAVIGATION 


This  section  describes  a  novel  fusion  and  learning  approach  developed  for  autonomous 
navigation  of  UAV  platforms.  The  approach  consists  of  incorporating  adaptation  into  a 
vector  integration  neural  navigation  model  in  conjunction  with  parametric  fuzzy 
adaptation  for  mission  fusion. 

2.1  Autonomous  Control  Methodology 

The  vector-integration-to-endpoint  model  of  the  neural  dynamics  of  planned  arm 
movements  consists  of  two  coupled  equations  which  may  be  directly  applied  to 
autonomous  platfonn  navigation  and  formation  control.  These  equations  form  a  model 
basis  for  sequentially  updating  the  present  platform  coordinate  position  in  real-time  to 
attain  a  desired  target  coordinate  for  platform  motion.  These  position  updates  are  driven 
by  a  mismatch  between  the  present  coordinate  position  of  the  platform  and  the  desired 
final  coordinate  target  position.  In  its  basic  form,  the  VITE  model  is  a  real-time  neural 
network  which,  for  the  \'h  entry  of  the  coordinate  vector,  may  be  expressed  as  [17],  [14]: 


=a(-F1+7;-/>) 

at 


(i) 


^pi:=G[  max(T„0)] 
at 


(2) 
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where  T,  is  the  i';'  entry  of  the  desired  position  coordinate  vector  T  of  the  platform;  Pt  the 
\lh  coordinate  of  the  present  position  command  coordinate  vector,  P\  Vt  the  \lb  coordinate 
of  the  mismatch  vector,  V;  and  G  is  a  possibly  time  varying  velocity  command  termed  the 
"GO"  signal. 

The  first  equation  generates  a  nonzero  mismatch  vector  V  whenever  the  present 
position  command  coordinate  vector  P  differs  from  the  desired  position  vector  T.  The 
mismatch  vector  is  continually  added  to  the  present  position  command  vector,  P,  until  the 
present  position  command  vector  equals  the  desired  position  vector  and  a  zero  valued 
mismatch  vector  results. 

The  above  trajectory  generator  permits  multiple  platforms  to  operate 
synchronously  in  time  even  though  the  coordinate  distances  from  present  position  to  final 
position  may  be  different  for  each  coordinate.  Furthermore,  the  gain  G,  the  "GO”  signal, 
regulates  the  speed  at  which  the  path  to  the  desired  target  position  evolves.  Even  though 
G  varies  with  time  synchronous  platform  operation  is  maintained.  Thus,  the  VITE  model 
is  capable  of  exhibiting  many  types  of  arm  motion  characteristics.  These  include  the 
capacity  to  speed  up,  slow  down  and  fast  freeze  anywhere  along  a  trajectory  path.  This  is 
accomplished  by  simply  adjusting  the  GO  signal  gains  in  the  model  since  the  value  of 
these  gains  determines  the  rate  at  which  the  trajectory  evolves  in  real-time.  Increasing 
the  gains  results  in  faster  motion,  decreasing  them  slows  motion  and  setting  them  to  zero 
causes  trajectory  evolution  to  freeze.  Such  a  feature  is  attractive  in  an  environment  where 
obstacles  may  suddenly  appear  along  a  trajectory. 

The  gain  G  also  provides  a  means  of  regulating  the  update  rate  of  the  present 
position  vector,  P,  by  employing  feedback  from  position  sensors  in  the  platform.  The 
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gain  G  is  adjusted  to  maintain  a  desired  tracking  error.  If  the  actual  tracking  error 
exceeds  that  desired,  the  gain  G  may  be  reduced  thereby  sacrificing  speed  and  improving 
accuracy.  On  the  other  hand,  if  the  tracking  error  is  less  than  that  desired,  the  gain  G  may 
be  made  larger  to  increase  speed  while  sacrificing  accuracy. 

In  order  to  assure  satisfactory  UAV  navigation  path  tracking  accuracy,  the 
discrepancy  between  trajectory  generator  present  position  commands  and  the  actual 
position  of  the  autonomous  platform  must  be  established  and  employed  by  the  trajectory 
generator.  A  structure  for  an  automatic  learning  process  to  accomplish  this  employs  error 
signals  derived  from  discrepancies  between  the  autonomous  platform  actual  position 
vector  and  that  seen  by  the  trajectory  generator.  This  provides  a  means  to  adaptively  tune 
the  present  position  vector  until  it  coincides  with  the  actual  position  vector  and  the 
trajectory  generator  learns  the  actual  position  of  the  platform.  This  type  of  learning 
happens  during  periods  of  tracking  inaccuracy,  i.e.  mismatch  and  allows  the  navigation 
updates  to  anticipate  platform  navigational  dynamics.  The  learning  and  trajectory  update 
functions  may  be  characterized  by: 

4^  =G[max(F/  0)]  +  G,[m*x(M„0)]  (3) 

at 

(4) 

at 

—  zi  -  SGp  {-£Zj  +  ma x(Mi  ,0)}  (5) 

dt 
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Equation  (3)  has  replaced  (1)  and  the  second  term  on  the  right  represents  a  gate  update 
signal  which  is  turned  on  by  the  gain  Gp.  Mx  is  another  mismatch  term  driven  by  any 
difference  between  the  actual  position,  yx,  and  the  present  position  command,  P,.  The 
variable  z;  is  an  associative  weight  which  readjusts  the  trajectory  generator  present 
position  until  it  coincides  with  the  actual  position.  The  gains  G  and  Gp  enable  external 
intelligent  control  of  the  update  and  learning  functions. 


2.2  Sensor  fusion  model  based  on  Navigation  control  methodology 

In  automating  the  operation  of  a  UAV  and  cluster  of  UAV’s,  a  fundamental  objective  is 
to  emulate  the  action  of  a  skilled  human  operator.  The  sensor  complement,  in  effect,  must 
allow  the  autonomous  intelligent  navigation  controller  to  discern  the  condition  of  the 
UAV  platform  system  and  to  a  degree  of  resolution  comparable  to  that  of  the  human 
operator.  This  includes  the  estimation  of  location  and  magnitude  of  piping  leakage 
failures,  assessment  of  sensor  functionality  and  failure  contingency  operation.  This 
involves  the  fusion  of  sensor  information  subject  to  constraints  imposed  by  the  system. 
This  may  be  viewed  as  a  constrained  estimation  and  identification  problem  which 
projects  noisy  sensor  data  onto  a  constraint  surface  representing  the  physical  model  of 
pressure  wave  propagation  [13],  [12],  [16],  The  intelligent  controller  assesses  the  action 
that  must  be  taken  to  optimize  system  performance.  This  requires  the  integration  of  the 
sensor  fusion  process  into  the  control  structure.  Two  methods  of  accomplishing  this 
integration  process  are  fuzzy  logic  and  neural  networks.  These  integration  approaches  are 
now  discussed. 
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2.3  Fuzzy  Theoretic  formulation  of  Automation 

The  process  of  fuzzy  logic  provides  a  means  of  establishing  a  projection  from  a  domain 
space  of  system  parameters  to  a  range  space  of  assessment  and  control  variables.  The  set 
of  system  domain  space  parameters  may  be  quite  diverse  including  not  only  various 
physical  variables  (e.g.  navigational  pipes,  temperatures,  propagation  speeds,  and 
pressures)  but  may  incorporate  variables  which  indicate  the  urgency  with  which  a  desired 
command  should  be  executed  by  the  system.  The  space  of  all  domain  variables  forms  the 
universe  of  discourse  for  fuzzy  inference.  The  domain  space  is  projected  or  mapped  into 
the  range  space  by  interpreting  a  set  of  rules.  The  rules  are  linguistic  statements  having  an 
antecedent-consequent,  i.e.  if-then,  fonnat.  The  rules  are  derived  from  an  intelligent 
understanding  of  the  operation  of  a  particular  system  and  can  reflect  imprecision  and 
uncertainty.  The  rule  set  is  evaluated  in  parallel  and  the  order  of  evaluation  is  not  critical. 
The  resulting  evaluation  of  the  rule  base  is  a  subset  set  of  elements  in  the  range  space.  In 
the  present  context,  the  range  space  elements  represent  condition  and  control  variables 
which  are  to  guide  the  system  in  executing  a  prescribed  command  in  the  "best  way" 
possible  subject  to  the  current  conditions  [13],  [23],  [24]. 

Prior  to  implementing  the  fuzzy  inference  process  for  a  particular  system  the 
variables  of  the  domain  and  range  space  must  be  defined  and  a  means  of  characterizing 
the  values  these  variables  can  take  on  must  be  provided.  In  addition  a  structure  for  the 
rule  base  must  be  established.  These  items  and  concerns  are  addressed  next  in  Section 
2.3.1  through  2.3.7.  In  addition  summary  overheads  are  provided  in  APPENDIX  B. 
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2.3.1  Fuzzy  Sets 


Fuzzy  sets  are  the  domain  and  range  elements  on  which  the  rule  base  operates.  Consider  a 
classical  set  such  as  A  =  {x|0  <  x  <  6}.  The  membership  of  this  set  is  well  defined  and  for 

any  x  e  R  either  xe  A  or  x  g  A  .  The  notion  of  a  fuzzy  set  stems  from  relaxing  the  sharp 
classification  criterion.  Membership  in  a  fuzzy  set  is  characterized  by  specifying  an 
element  along  with  a  measure  of  how  much  the  element  belongs  to  the  set.  That  is,  the 
elements  of  a  fuzzy  set,  F ,  have  the  ordered  pair  form  (x,  m(x))  e  F  where  x  a  well- 
defined  element  as  in  the  classical  set  case  and  m(x)  is  the  membership  function 
associated  with  how  strongly  x  fits  into  F .  The  values  of  a  membership  function  are 
typically  restricted  to  the  range  0  <  m(x)  <1.0.  Where  m(x)  =  0  indicates  that  x  is  least 
likely  to  fit  into  the  fuzzy  set  while  m(x)  =  1.0  indicates  that  x  most  likely  to  be  a 
member  of  the  fuzzy  set.  This  set  structure  allows  two  things  to  happen.  Firstly  all 
subsets  of  the  domain  space,  which  correspond  to  rule  premises,  have  membership  of  an 
element  x  assessed  by  a  mapping  to  a  membership  value  between  0.0  and  1.0.  Secondly 
these  mapping  values  provide  a  means  of  projecting  how  strongly  the  premises  of  a  rule 
are  satisfied  and  hence  how  strongly  the  conclusion  should  be  emphasized. 

2.3.2  Fuzzy  Logical  Operations 

The  interrelationship  between  the  elements  and  subsets  in  classical  set  theory  are  well 
defined  by  an  algebraic  structure  terms  of  unions,  intersections,  complements  and 
differences.  These  operations  are  closed,  commutative,  associative  and  distributive.  The 
truth  or  validity  of  a  logical  statement  may  be  determined  from  these  operations  and 
properties.  An  extension  of  this  process  has  been  developed  for  fuzzy  set  algebra.  Let 
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A  =  {(x,m/ ,(x))}  and  B  =  {(x,  mB  (x))}  be  fuzzy  sets  over  the  domain  space.  Then  the 
union  of  these  sets  is  given  by  Ar\B  =  {x,S(mA(x),mB(x)} .  Where  the  properties  of  the 
function  S  are  chosen  to  yield  classical  set  theoretic  relationships  when  the  membership 
functions  are  reduce  the  binary  values  of  0  and  1 .  Thus  classical  sets  are  a  subspace  of 
fuzzy  sets.  The  function  S  must  be  monotonic,  commutative  and  associative.  Two 
commonly  employed  functions  which  satisfy  these  properties  are  the  minimum  function 
relation  given  by  S  =  min^n^  (x),mB(x))  and  the  product  function  relation  given  by 
S  =  product{mA(x),mg(x)).  Further  continuing  in  the  same  vein,  the  union  of  two  fuzzy 
domain  space  sets  is  given  by  A  'O  B  =  {x,  5'(w/)(x),  m6(x)}  where  S  again  must  satisfy 
the  conditions  stated  above.  Commonly  employed  functions  which  satisfy  these 
conditions  include  the  maximum  relation  given  by  S  =  max(/n/l(x),mjg(x))and  also  the 
relation  5  =  mA  (x)  +  mB  (x)  -  mA  (x)  x  mB  (x)  which  has  been  termed  the  probabilistic  OR 
function. 

2.3.3  Fuzzy  Logic  Process 

In  implementing  a  fuzzy  inference  system  several  processes  must  be  accomplished  to 
map  the  domain  space  variables  to  the  output  space.  These  processes  consist  of:  (1)  a 
fuzzification  interface  for  the  input  variable;  (2)  applying  fuzzy  union  and  intersection 
operations;  (3)  applying  a  fuzzy  inference  mechanism;  (4)  defuzzification  interface  which 
maps  the  conclusions  of  the  inference  mechanism  to  generate  an  output  with  the  objective 
of  controlling  and  or  assessing  the  system  [30],  [24],  [21].The  method  in  which  these  four 
steps  or  processes  are  implemented  characterized  in  large  part  the  response  dynamics  of 
the  fuzzy  deduction  process.  Further  details  about  each  of  these  steps  are  considered  next. 
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2.3.4  Fuzzification 


The  process  of  fuzzification  consists  of  assigning  the  degree  to  which  each  parameter  in 
the  domain  space  belongs  to  the  fuzzy  sets  which  comprise  the  premises  for  the  deduction 
rule  set.  This  determines  which  rules  are  to  be  used  and  the  extent  to  which  the  rule  will 
influence  the  final  outcome  of  the  fuzzy  logic  process.  Fuzzification  is  accomplished  by 
employing  a  membership  function.  The  fundamental  characteristics  of  a  sigmoidal 
membership  function  are  shown  Figure  2.  The  range  of  the  membership  function  is  from 
0.0  to  1 .0.  This  range  of  values  expresses  the  confidence  with  which  domain  space 
elements  belong  to  a  particular  fuzzy  set.  The  membership  function  value  1 .0  represents 
certain  membership  while  the  membership  function  value  0.0  indicates  lack  of 
membership.  The  process  of  Fuzzification  consists  of  resolving  the  input,  i.e.  the  domain, 
variables  into  their  degree  of  membership  in  each  of  the  antecedent  fuzzy  sets.  In  the  case 
where  the  fuzzy  rule  base  contains  rules  with  multiple  part  antecedents,  fuzzy  set 
operations  must  be  applied  to  blend  the  antecedents  into  a  single  number  between  0.0  and 
1.0.  This  process  typically  involves  the  intersection  (AND)  and  union  (OR)  operations  on 
fuzzy  set  elements  previously  defined.  In  any  case,  the  fuzzification  process  assigns  a 
single  number  to  the  antecedent  for  each  rule  in  the  rule  base.  Then  next  step  in  the 
process  is  one  of  assessing  the  validity  of  the  conclusions  associated  with  the  various 
rules. 
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Figure  2:  Illustration  of  a  membership  function  for  a  domain  space  variable  x. 

2.3.5  Rule  Implication  Process 

The  process  of  implication  determines  what  the  conclusion  of  the  fuzzy  logic  process  will 
be.  The  previous  fuzzification  process  determined  which  rules  are  active  in  the  inference 
process.  The  recommendation  of  each  rule  is  first  considered  individually.  Each  rule  has 
associated  with  it  a  weight  between  0.0  and  1.0  which  is  applied  to  the  antecedent.  This 
rule  weight  provides  a  means  of  tuning  the  relative  influence  of  each  rule  on  the  outcome. 
These  weighted  individual  recommendations  are  then  blended  into  an  aggregate 
conclusion  which  is  the  net  recommendation  of  the  fuzzy  inference  system  in  controlling 
the  system  or  estimating  its  condition. 

The  outcome  of  each  rule  is  a  consequent  that  is  a  fuzzy  set  with  an  associated 
membership  function  that  weighs  the  linguistic  content  of  the  rule.  The  input  to  the 
implication  process  is  a  single  number  derived  from  the  antecedent  and  the  resulting 
output  is  a  fuzzy  set.  The  implication  process  is  done  on  a  premise  fuzzy  set  derived  by 
an  intersection,  i.e.  an  AND,  operation  where  the  respective  membership  functions  are 
derived  from  5  s  min  ,  i.e.  minimum,  operator  or  from  S  =  prod  ,  i.e.  product,  operator. 
The  former  operation  truncates  the  resultant  fuzzy  set  while  the  latter  operation  scales  the 
resultant  fuzzy  set. 
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2.3.6  Aggregation  of  the  Rule  Outputs 

The  fuzzy  inference  system  in  effect  tests  all  of  the  rules  in  the  rule  base  and  then  must 
combine  all  of  the  rule  results  in  some  fashion  so  as  to  derive  a  net  decision  or 
conclusion.  The  process  of  aggregation  consists  of  combining  the  fuzzy  sets  that  are  the 
outputs  of  each  rule  in  the  rule  base  to  form  a  single  fuzzy  set.  The  resulting  single  fuzzy 
set  represents  the  final  conclusion  of  the  entire  rule  base.  The  aggregation  method  is 
typically  a  commutative  process  in  which  the  order  in  which  the  rules  are  applied  is 
unimportant.  The  fuzzy  sets,  which  are  the  consequence  of  the  rule  base,  may  be 
combined  commutatively  with  several  common  functions.  The  membership  functions  of 
the  rule  consequents  may  be  combined  using  the  maximum  function,  the  probabilistic 
OR,  and  the  sum  function.  The  aggregation  process  employed  in  this  work  is  the 
maximum  function  as  described  in  the  previous  section.  Application  of  this  function 
commutatively  combines  all  of  the  membership  functions  from  the  consequent  of  each 
rule  into  an  aggregate  membership  function  that  describes  the  resulting  final  fuzzy  set, 
i.e.  the  conclusion  derived  from  the  entire  rule  base. 

2.3.7  Defuzzification  to  Generate  Output 

The  structure  of  the  controller  considered  requires  that  the  output  of  the  fuzzy  logic 
process  be  a  scalar  quantity  as  opposed  to  a  fuzzy  set.  Defuzzification  is  the  conversion 
of  a  fuzzy  set  characterization  of  the  outcome  of  the  fuzzy  logic  process  to  a  single  scalar 
quantity.  The  output  of  a  fuzzy  logic  process  consists  of  the  logical  union  of  two  or  more 
fuzzy  membership  functions  defined  on  the  output  variables.  The  center-of-gravity 
method  is  employed  to  defuzzify  the  output  of  the  fuzzy  logic  process.  This  method  is 
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perhaps  the  most  well  known  and  physically  tractable  approach  to  defuzzification,  [24]. 
A  scalar  output,  unel,  is  obtained  by  operating  on  the  membership  function,  mllel(x) ,  of 
the  conclusion  of  the  fuzzy  logic  process.  The  scalar  output  is  given  by  the  expression 


(6) 


u, 


net 


where  x  represents  the  range  of  values  which  the  fuzzy  set  output  encompasses.  Thus  the 
center-of-gravity  method  calculates  an  average  unel  of  all  possible  output  values  x 

weighted  by  the  degree  of  membership  mnel{x) of  the  output  fuzzy  set.  This  method  is 

heuristic  and  generally  performs  well  if  the  experience  of  system  operators  is  integrated 
with  analytical  system  models. 

2.4  Neural  network  sensor  fusion  and  autonomous  control  action 

An  architecture  for  UAV  system  intelligent  autonomous  control  employing  a  neural 
network  is  shown  in  Figure  3.  One  of  the  most  crucial  considerations  in  formulating  a 
neural  based  control  process  is  the  selection  of  neural  network  inputs  and  outputs  [14], 
[15],  [21],  [30],  In  the  proposed  neural  controller,  navigation  position  readings  along  with 
navigation  generator  states  have  been  selected  as  inputs  to  the  network.  The  network 
outputs  are  the  various  steering  and  thrust  control  signals.  The  neural  network  has  the 
capability  to  leam  dynamic  control  sequences  and  to  interpolate  intelligent  control 
actions  for  system  states  which  were  not  part  of  the  original  training  algorithm. 
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Figure  3:  A  neural  network  architecture  for  autonomous  UAV  system  control. 

Another  crucial  factor  in  the  implementation  of  the  neural  network  is  the  training 
sequence  selected.  The  neural  network  can  provide  significant  advantages  in  cases  where 
the  system  model  is  not  well  defined  and/or  where  a  system  inverse  controller  is 
employed.  In  these  cases  the  network  may  be  trained  from  data  generated  from  actual 
system  operation.  In  cases  where  the  analytical  model  of  the  system  is  well  defined,  the 
neural  network  may  provide  a  means  of  realizing  various  control  actions  which  may  be 
non-linear.  In  the  context  of  the  UAV  system  an  autonomous  intelligent  control  rule  and 
action  base  must  be  established  prior  attempting  a  training  algorithm  formulation. 

2.5  SIMULATION  MODEL 

The  baseline  considerations  in  using  autonomous  platforms  are:  navigational  control  of 
autonomous  platforms  and  platform  clusters;  sensing  capability  for  environment/mission; 
and  the  intelligent  fusion  of  sensor,  mission  and  control  data.  Key  to  the  success  of  this 
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concept  is  the  ability  to  maintain  functionality  while  restricting  navigational  and  sensor 
information  bandwidths  to  acceptable  levels.  A  vector  integration  to  endpoint  navigation 
model  is  employed  as  a  distributed  automata  agent  to  satisfy  the  operational  and 
information  rate  concerns.  Simulation  results  illustrate  the  different  characteristics  of  this 
model  based  on  equations  (1)  and  (2). 

2.5.1  Autonomous  Vehicle  Challenges 

The  main  challenges  in  implementing  a  navigational  system  include  realizing  the 
properties  of: 

1)  Fast  real-time  capability: 

a)  Maintain  computational  efficiency  of  the  navigational  model. 

b)  Minimize  communication  bandwidth  between  different  units. 

2)  Implicitly  provide  for  the  fusion  of  mission  objectives  and  navigation  goals.  The 
goal  of  the  above  properties  is  the  achievement  of  tracking  accuracy  as  which  seeks  to 
minimize  the  difference  between  trajectory  generator  present  position  and  the  actual 
position  of  the  autonomous  platform.  The  rate  of  navigation  progression  must  be 
adjusted  based  on  the  ability  of  the  autonomous  platform  to  follow  trajectory  commands. 
In  effect,  the  navigation  engine  learns  the  UAV  response  characteristics  and  adjusts 
navigation  commands  to  accommodate  them. 

2.5.2  Obstacle  Avoidance  Rule  Set 

One  of  the  primary  navigational  functions  which  must  be  accomplished  is  that  of  obstacle 
avoidance.  The  navigation  model  must  steer  the  vehicle  to  clear  any  obstacles  with 
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respect  to  an  x,  y  coordinate  system,  the  baseline  rule  set  chosen  for  collision  avoidance 
is  as  follows: 

-  If  Ax  is  large  then  make  no  correction. 

-  If  Ay  is  large  then  make  no  correction. 

-  If  Ax  is  small  and  Ay  is  small  then  take  corrective  action. 

-  If  Ax  is  small  and  Ay  is  smallest  then  take  corrective  action  emphasizing  y. 

-  If  Ax  is  smallest  and  Ay  is  small  then  take  corrective  action  emphasizing  x. 

2.5.3  Navigation  Model  Simulation  Performance 

The  performance  of  the  integration  to  endpoint  based  navigation  model  is  presented  for 
the  cases  of  open  sea  navigation  with:  (1)  no  hazards;  (2)  with  hazards;  (3)  with  moving 
hazards  and  (4)  finally  intercept  of  moving  targets  without  and  with  stand  off. 

2.5. 3.1  Open  Sea  Navigation  Without  Obstacle 

Figure  4  shows  the  simulated  path  motion  of  the  vehicle  with  no  obstacle  present  in  its 
path.  The  vehicle  starts  from  coordinates  (2,2)  with  its  final  target  position  being  (10,10). 
The  navigation  motion  is  inherently  smooth  in  position  and  velocity.  Also,  the  changes  in 
x  and  y  coordinates  are  synergetic  from  start  to  stop.  The  generator  exhibits  non- 
oscillatory  convergence  throughout  the  open  sea  field  of  navigation. 
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Figure  4:  Navigation  path  without  obstacles 

2. 5. 3. 2  Open  Sea  Navigation  With  Hazards 

Figures  5  and  6  illustrate  typical  simulation  results  for  the  motion  of  a  vehicle  around 
circular  hazards  of  differing  diameters.  The  vehicle  autonomously  changes  its  path  as  it 
nears  obstacles.  It  is  observed  that  the  coordinate  points  become  more  dense  near  the 
final  position  indicating  smooth  slowing  velocity  as  the  coordinates  are  approached.  The 
navigation  trajectories  generated  continue  to  exhibit  non-oscillatory  monotonic 
convergence  which  is  smooth  in  position  and  velocity  throughout  the  navigation  region. 
A  functional  characterization  of  the  obstacle  avoidance  navigation  engine  is  presented  in 
Figure  7. 


18 


Figure  5:  Navigation  path  with  obstacle  at  (4,4) 


Figure  6:  Navigation  path  with  obstacle  of  diameter  1.16  at  (4,4) 
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Figure  7:  Functional  characterization  of  obstacle  avoidance 
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2. 5.3.3  Open  Sea  Navigation  with  Moving  and  Fixed  Hazards 

Figures  8  and  9  show  two  vehicles  on  a  collision  course.  In  this  case  the  obstacle 
avoidance  algorithm  has  been  modified  to  accept  moving  obstacle  with  each  vehicle 
constituting  a  moving  obstacle  for  the  other  vehicle.  As  can  be  seen  in  Figure  8,  one  of 
the  vehicles  autonomously  and  independently  slows  down  so  as  to  avoid  UAV  to  UAV 
collision. 

Figure  9  illustrates  a  navigation  scenario  where  both  moving  and  fixed  navigation 
hazards  are  present.  Here  the  navigation  generator  autonomously  acts  to  simultaneously 
avoid  both  the  fixed  and  moving  obstacles. 


Figure  8:  Navigational  path  with  a  moving  obstacle 

The  generated  trajectories  are  again  smooth  in  position  and  velocity  exhibiting  monotonic 
synergistic  convergence  to  target  coordinates. 
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Figure  9:  Navigational  path  with  a  fixed  and  moving  obstacle 

2. 5.3.4  Pursuit  and  Interception 

The  interception  of  a  moving  target  is  a  common  navigation  maneuver.  This  scenario  is 
shown  in  Figure  10.  In  addition,  if  the  intercepted  target  has  sensing  or  defense  capability 
it  may  be  necessary  to  follow  the  target  platform  with  a  stand  off  distance. 

This  figure  demonstrates  inherent  capability  of  the  navigation  generator  to  intercept  a 
moving  target  while  maintaining  a  safe  stand-off  distance.  The  simulation  results  of 
stand-off  pursuit  are  given  in  Figure  1 1 . 
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Figure  10:  Moving  target  intercept  navigational  path 


Figure  11:  Moving  target  intercept  with  standoff 


2.5.4  MIMO  Real  Time  Model 


A  MIMO  model  was  derived  from  Equations  (1)  and  (2)  for  the  purpose  of  allowing 
navigational  computation  to  be  interleaved  with  real-time  sensing  and  communication.  It 
also  provides  a  means  of  fusing  mission  and  navigation  information  as  updated  by 
communications  information.  In  addition,  a  MIMO  Real-time  navigation  model  allows 
for  the  sequential  processing  of  navigation,  GPS  and  communication  data  in  the 
computation  of  a  neuro-fuzzy  PIM.  The  baseline  real  time  implementation  of  the  MIMO 
sequential  mode  was  accomplished  first  with  communication  links  consisting  of  serial 
cables  and  simulated  GPS  position  words.  This  scenario  is  described  in  the  next  section. 

2.5.4. 1  Real  Time  Performance  &  Emulation 

The  real  time  setup  includes  three  computers;  two  slave  computers  and  one  master 


RS  232  Cables 


Slave  1  Master  Slave  2 


Figure  12:  Set  up  for  real  time  simulations 

computer.  The  units  are  all  connected  by  means  of  RS  232  cables.  In  the  real-time  setup 
both  the  slave  computers  wait  for  mission  start  up  in  addition  to  mission  commands  sent 
as  targets  and  obstacles.  The  master  station  controls  and  monitors  the  mission  progress. 
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2.5.4. 1.1  Real  Time  Collision  Avoidance 


Once  the  mission  commences  Real-time  collision  avoidance  is  characterized  by 
live  communication  link  which,  assumes  that  both  slave  computers  provide  their  position 
coordinates  to  the  master  computer  as  well  as  to  the  other  slave  computer.  The  two  slave 
computers  autonomously  navigate  their  respective  UAV’s  in  response  to  the  location  of 
the  other  UAV  and  the  coordinates  of  mission  targets.  Figure  13  illustrates  typical 
autonomous  navigation  results  obtained  during  real-time  navigation  emulation  for  the 
typical  case  of  UAV’s  with  two  obstacles  in  the  direct  path.  In  this  case  the  desired  final 
positions  are  communicated  from  the  master  station.  The  vehicle  is  observed  to  smoothly 
change  its  path  depending  on  the  position  of  the  obstacles. 


Figure  13:  Real  time  navigation  path  with  two  obstacles 
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A  real-time  navigation  case  with  two  moving  vehicles  is  presented  next.  In  the 
real-time  set  up,  motion  of  both  the  vehicles  can  be  seen  on  the  same  plot  at  the  master 
station.  Figure  14  presents  the  real-time  navigation  emulation  results.  Both  vehicles  are 
observed  to  smoothly  avoid  both  moving  and  fixed  obstacles.  All  mission  and  obstacle 
data  are  communicated  via  serial  ports  and  the  resulting  paths  generated  are  non- 
oscillatory  and  monotonically  convergent. 
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Figure  14:  Plot  of  a  real  time  navigational  path  with  a  fixed  and  moving  obstacle 

(screenshot  at  master  computer) 

2. 5.4. 2  Real-Time  Performance 

The  observed  real-time  performance  indicates  that  the  trajectory  paths  exhibit  a  fast  ramp 
at  start-up  and  a  slow  down  as  the  final  target  is  approached.  The  trajectories  are  stable, 
non-oscillatory,  monotonic  and  smooth  in  position  and  velocity.  Some  nonsymmetry  in 


26 


the  smoothness  was  observed  with  accelerated  initial  start  up  and  deaccelerated  final 
approach  to  target. 


2.6  Real  Time  Implementation  Approach 

The  considerations  in  implementing  the  previously  described  real-time  navigation 
generator  with  micro-controller  and  intelligent  control  is  now  addressed.  From  a  broad 
perspective,  the  development  of  innovative  technology  aimed  at  reducing  manning 
requirements  and  life  cycle  costs  are  driven  by  the  Smart/Intelligent  Ship  philosophy 
evolved  and  adopted  by  the  Navy  for  future  platforms.  The  development  of  autonomous 
navigation  control  and  information  integration  systems  is  anticipated  to  play  a  major  role 
in  this  initiative.  Within  the  framework  of  this  initiative,  the  development  of  a 
methodology  for  the  design  and  implementation  of  autonomous  machinery  control 
systems  is  proposed.  The  methodology  must  employ  available  COTS  technology  to 
provide  a  flexible  integrated  approach  to  accomplish  the  essential  tasks  of  modeling, 
design,  analysis,  and  assessment  for  shipboard  environment/mission,  prototyping  and 
real-time  implementation.  An  innovative  feature  of  the  proposed  approach  is  that  the 
previous  task  sequence  is  completely  integrated  into  the  final  design  implementation.  The 
result  is  that  all  developmental  expertise  is  retained  within  the  controller  and  readily 
available  for  future  use.  The  software  structure  also  allows  on-the-fly  contingency 
assessment,  parameter  adjustment  and  supervisory  operator  training. 
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2.6.1  Considerations  for  Autonomous  Machinery  Control 

The  development  of  a  flexible  methodology  for  the  integrated  design  and  implementation 
of  autonomous  intelligent  machinery  control  systems  is  proposed.  The  process  of 
developing  and  deploying  an  intelligent  autonomous  navigation  control  system  includes 
the  following  interrelated  phases:  (1)  formulation  of  control  objectives  and  strategy 
which  incorporates  component,  system  and  shipboard  wide  intelligence  capability;  (2) 
modeling;  (3)  control  algorithm  analysis,  design  and  assessment;  (4)  simulation;  (5)  rapid 
controller  prototyping  test  and  evaluation;  (6)  final  implementation  (distributed  in  COTS 
DSP  based  hardware  with  supervisory  information  integration).  Considerations 
associated  with  each  of  these  phases  are  discussed  below. 

2.6. 1.1  Formulation  of  control  objectives  and  strategy 

The  formulation  of  control  objectives  is  driven  by  the  need  for  any  autonomously 
operated  system  to  reconfigure  itself  to  best  support  the  mission  of  the  platform.  This 
includes  anticipated  modes  of  normal  operation  and  damage  control  subject  to  overall 
ship  status  and  immediate  mission  needs.  In  addition,  the  development  of  an  intelligent 
decision  structure  is  desirable  in  an  attempt  to  provide  reconfiguration  capability  under 
operational  conditions  resulting  from  unanticipated  scenarios. 

2. 6. 1.2  Modeling 

Modeling  is  the  essential  fundamental  process  upon  which  control  system  analysis  and 
design  is  based.  Time  and  frequency  domain  models  are  based  on  an  understanding  of  the 
physical  processes  that  characterize  the  system  components.  These  yield  analytical  and 
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numerical  based  approximate  quantification  of  system  component  response.  Component 
models  may  include  the  plant,  associated  sensors  and  actuator  dynamics. 

2. 6. 1.3  Simulation 

Simulation  provides  the  essential  vehicle  to  integrate  component  models  into  a  complete 
system  model  that  also  includes  the  controller.  Closed  loop  simulations  allow  an 
assessment  of  how  well  the  initial  controller  structure,  and  subsequent  iterations  of 
controller  design,  meet  performance  specifications. 

2.6. 1.4  Distributive  controller  design,  assessment  and  algorithm 

The  framework  for  controller  development  is  one  of  autonomous  distributive  control 
which  is  capable  of  adapting  to  numerous  shipboard  operational  scenarios,  anticipated 
and  unanticipated,  realizing  a  fight-whil e-hurt  capability.  The  distributive  control  process 
incorporates:  micro-controller  component  level  intelligence;  and  system  level  intelligence 
via  corroborative  micro-controller  components  and  PC  based  software.  In  addition,  the 
control  architecture  must  incorporate  supervisory  human  direction  with  the  potential  to 
accept  threat  status  and  evasion  information  from  combat  systems,  navigation, 
communication  and  damage  control  shipboard  functions.  The  process  of  controller  design 
is  twofold.  The  first  component  is  control  strategy  development  and  commonly  involves 
a  blend  of  approaches  which  include  classical,  modem,  robust,  p-control, 
adaptive/predictive,  intelligent,  optimal,  neural  and  fuzzy  methods.  The  second 
component  deals  with  implementation  issues.  Once  a  control  algorithm  is  specified  it 
must  be  encoded  in  a  distributive  environment  consisting  of  series  DSP  based  micro- 
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controllers.  Consideration  must  be  given  to  processor  execution  speed,  I/O  buffering, 
A/D  and  D/A  conversion,  and  synchronization.  From  the  distributive  perspective  issues 
of  concern  in  best  satisfying  shipboard  platform  mission  needs  include  the  stand  alone 
control  action  of  individual  micro-controllers  in  concert  with  corroborative  control 
actions  of  micro-controller  subsets. 

2. 6.1. 5  Control  algorithm  prototype  implementation  and  real-time  testing 
No  matter  how  refined  the  model  of  an  autonomous  system  is,  it  represents  the  result  of 
an  abstraction  process  and  consequently  there  are  physical  phenomena  that  are  not  totally 
accounted  for  by  the  model.  The  model  serves  as  the  vehicle  for  bringing  to  bear  the 
theoretically  powerful  analytical  analysis  and  design  approaches  to  the  problem  of 
developing  an  optimized  automated  machinery  system.  The  resulting  design  is  optimized 
with  respect  to  the  model  representing  the  system  and  its  validity  in  the  framework  of  an 
actual  hardware  implementation  must  be  addressed. 

The  conclusive  performance  assessment  of  a  control  algorithm  lies  in  real-time 
implementation  testing.  This  configuration  consists  of  a  physical  plant,  full  or  scaled,  and 
in-the-loop  controller  hardware  in  addition  to  the  implementation  of  supervisory 
software/hardware  to  assess  and  tune  the  control  algorithm  in  real  time.  The  supervisory 
software/hardware  may  also  serve  as  a  means  to  evaluate  and  optimize  available  control 
options  subject  to  condition  assessment,  mission  and  human  operator  inputs. 
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3.0  DISTRIBUTRED  DIVERSE  SENSING 


In  the  conventional  array  signal  processing  approach,  the  sensor  locations  must  be 
separated  by  a  maximum  of  one  half  wavelength  of  the  dominant  signal  component  in 
order  to  estimate  the  direction  of  arrival  of  acoustic  source  waves.  This  scenario  is 
applicable  to  a  sensor  node  containing  a  micro-array  of  acoustic  sensors.  In  the  case 
where  the  UAV  sensor  nodes  contain  single  acoustic  sensors  and  are  separated  by  more 
than  a  half  wavelength  of  the  dominant  acoustic  frequency  component,  an  augmented 
sensor  processing  problem  results.  The  augmented  sensor  problem  is  illustrated  in  Figure 
15.  In  this  case  the  signal  model  cannot  be  directly  employed  to  deduce  source  location 
information  due  to  spatial  sparsity.  This  may  be  overcome  however,  by  employing  the 
platform  model  dynamics,  which  augment  the  signal  samples  with  a  physical  model  of 
the  propagation  and  acoustic  environment.  The  augmented  modeling  cases,  which  were 
considered,  include  the  sensor-to-sensor  feature  propagation  delay  augmentation  and  the 
hyper-geometric  feature  propagation  delay  augmentation.  Geolocation  and  bearing 
technologies  currently  are  rapidly  evolving  and  the  potential  exists  for  employing  these  in 
extending  the  augmented  location  problem  to  a  more  autonomous  state. 
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Figure  15:  Framework  of  the  augmented  sensor  array  problem. 

3.1  Sensor-to-sensor  feature  propagation  delay 

Propagation  delay  source  location  is  characterized  by  sensor  distributions  that  may  be 
randomly  distributed  throughout  a  region  without  any  restrictions  on  the  distances 
between  sensors  in  relation  to  the  acoustic  source  wavelength  [19],  [26].  This  is  in 
contrast  to  phase  based  conventional  sensor  array  processing  techniques,  which  bound  the 
spacing  between  sensors  to  within  a  half  wavelength  of  the  dominant  frequency  being 
analyzed  for  direction  of  arrival  estimation.  The  typical  geometric  layout  of  the  feature 
propagation  approach  is  illustrated  in  Figure  16.  The  respective  propagation  delay  times 
are  given  by: 
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(7) 


=  Ti~T\ 

at{J=t3-t, 

ATU  =Ta-Tx 

where  7J ,  T2 ,  T3  and  7^  are  the  times  of  arrival  at  the  corresponding  sensors.  The  relative 

propagation  delay  between  sensors  of  the  source  feature  is  used  in  concert  with  the 
coordinates  of  the  acoustic  sensor  elements  to  determine  the  location  of  the  acoustic 
source  within  the  sensing  region.  Considering  the  source  to  have  coordinates  u  and  v,  the 
equations  of  concentric  circles  about  the  source  with  radii  of  lengths  corresponding  to  the 
various  distances  from  the  source  to  a  respective  sensor  may  be  expressed  as 
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(8) 


(x2  -uf  +0,  -v)2  =(fl,)2 

( x2  ~u )2  +(j;2 -v)2  =  (^i  +^Ar,2)2 
(x3-W)2+(j3-v)2  =(R,+sATu)2 
(x4-u)2+(y4-v)2=(Rl+sAT[4)2 


Equation  (8)  may  be  rearranged  into  the  matrix  form  as 


2xj 

-2*2 

2  T. 

2y  2 

-2vA7’12 

u 

v2A7’12  +x2  +  y2  —x\  —  _y2 

2xj 

-2.v3 

2f. 

-2y3 

—  2vA7j3 

V 

= 

v2A  Tx\+x2+y2-xl-yl 

2xj 

-2x4 

2y, 

-2  y4 

-2vA  rI4_ 

_v2AT2+x2+y2-xl-y24_ 

The  matrix  equation  (9)  may  be  directly  solved  to  yield 


u 

V 


*1 


and  requires  a  minimum  of 


four  acoustic  sensors.  This  method  is  may  be  extrapolated  to  three  dimensions  and 
requires  a  minimum  of  five  non-collinear  acoustic  sensors  to  achieve  an  estimate  of  the 
source  coordinates. 


3.2  Hyper-geometric  feature  propagation 

The  hyper-geometric  characterization  of  source  location  with  respect  to  a  constellation  of 
acoustic  sensors  is  illustrated  in  Figure  17.  In  contrast  to  the  four  sensors  required  by  the 
feature  time  delay  approach  for  2-D  source  location  estimation,  the  hyper-geometric 
approach  requires  only  three  sensors,  which  are  non-collinear.  In  the  case  of  3-D  source 
coordinate  estimation,  the  hyper-geometric  approach  requires  four  non-collinear  acoustic 
sensor  vice  the  five  non-collinear  sensors  required  by  the  feature  propagation  delay 
approach.  The  hyper-geometric  approach  consists  of  simultaneously  solving  for  the 
intersection  of  two,  2-D  case,  or  three,  3-D  case,  hyper  curves.  The  appropriate  curves  are 
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selected  based  on  the  time  of  arrival  of  the  acoustic  wave  at  the  various  sensors.  In  the  2- 


D  case  the  two  sensors  with  the  first  two  arrival  times  are  chosen  as  the  corresponding 
foci  of  hyper-geometric  curves.  In  the  three  dimensional  case  the  three  sensors  with  the 
first  three  arrival  times  are  chosen  as  the  reference  foci.  Without  loss  of  generality,  the 
hyper-geometric  method  may  be  characterized  for  the  case  of  sensors  symmetrically 
located  about  the  origin.  Let  the  lattice  rectum  intersect  the  x-axis  at  c  and  — c,  which  are 
foci.  These  curves  are  characterized  by 


x 


2 


a 


2 


9 


(10) 


In  the  analogous  case  of  focii  on  the  y-axis  at  c  and  -c,  the  corresponding  curves  are 
represented  by 


y 


2 


a 
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(ii) 
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Sensor  1 


Figure  17:  Hyper-geometric  feature  propagation  delay  sensor 

geometry. 


The  quantities  ay  and  ax  are  they  and  x-axis  intercepts  as  derived  from  a  fusion  of 
the  relative  sensor-to-sensor  time  delays.  Furthermore,  the  intersection  between  these  sets 
of  curves  may  be  expressed  as 
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(12) 


and 


where 


x2  = 


fe  +al) 

fbl  «.?' 


y 


y  =  - 


*L 


(13) 


(14) 


and 


(15) 


The  generalized  coordinate  case  is  derived  from  the  above  symmetric 
development  by  incorporating  applicable  rotational  and  transnational  operators  on  the 
corresponding  2-D  or  3-D  coordinate  n-tuple.  The  localization  process  requires  only  the 
extraction  of  the  applicable  ax,  aY  or  a,  from  the  sampled  acoustic  data.  Thus  yielding  a 
reduction  in  the  required  number  of  sensors  required  in  contrast  to  the  direct  feature 
propagation  delay  methodology. 
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3.3  DISTRIBUTED  SHIPBOARD  APPLICATION 


The  fundamental  problem  of  autonomous  intelligent  control  for  distributed  fluid  systems 
is  addressed  as  an  initial  target  application.  Modeling,  simulation  and  control 
methodologies  are  considered  for  leak  damage  detection,  leak  isolation  and  resource 
optimization  in  meeting  load  demand.  While  the  results  are  primarily  applied  to  UAV 
distributed  systems,  they  extend  directly  to  a  wide  range  of  sensors  and  deployment 
sensors. 

3.3.1  System  Configuration 

A  representative  low-pressure  (LP)  system  configuration  is  illustrated  in  Figure  18.  The 


_  Piping 

Figure  18: 


Micro-controller  with  sensors 


'J' 


Leakage  loss 


A  typical  UAV  system  series-parallel  network  configuration. 
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topology  of  the  system  consists  of  a  series-parallel  fluid  flow  network  composed  of 
piping  sections,  valves,  sensors,  micro-controllers  and  fluid  supplies.  The  fluid  supply 
units  feed  the  network  and  consist  of  a  motor-driven  regulated  compressor  unit,  a  filter 
unit  and  a  dehydration  baffle.  The  piping  sections  contain  automated  on/off  flow  valves 
whose  control  actions  are  dictated  by  a  dedicated  micro-controller  unit.  The  micro¬ 
controller  units  each  employ  a  locally  programmable  control  algorithm  driven  by  the 
readings  of  two  pressure  switch  sensors  located  on  either  side  of  the  valve  unit.  Thus,  a 
micro-controller  operates  its  valve  unit  in  accordance  with  an  assessment  of  the  pressure 
state  of  system  piping  on  both  sides  of  the  valve.  In  addition,  each  of  the  piping  sections 
and/or  valve  units  may  experience  a  leakage  failure.  Under  normal  and  failure  operating 
conditions,  the  fundamental  goal  is  to  achieve  optimal  autonomous  operation  in  meeting 
the  system  load  demand. 

3.3.2  Distributed  Failure  Problem  Statement 

The  problem  focus  is  one  of  developing  a  distributed  control  methodology  with  the 
objective  of  achieving  optimized  autonomous  operation.  The  major  aspect  of  this 
problem  is  how  to  allocate  the  available  capability  of  the  UAV  supply  to  best  satisfy  load 
demands.  This  includes  supplying  load  demands  under  normal  operation  and  subject  to 
anticipated  and  unanticipated  system  damage  scenarios.  The  system  component  damage 
may  include  any  combination  of:  (1)  moderate  leakage  failures  in  piping  and  valve 
components;  (2)  catastrophic  leakage  failures  in  piping  and  valve  components;  (3)  stuck 
valve  failure;  (4)  partial/total  failure  of  communication  between  micro-controller  units  (5) 
gradual  component  deterioration. 
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3.3.3  Approach 

The  approach  taken  in  addressing  the  stated  problem  consists  of:  (1)  a  characterization  of 
the  dominant  physical  phenomena  within  the  system;  (2)  analytical  component  model 
development;  (3)  simulation  component  model  development  in  conjunction  with  micro 
controller  action  modeling;  (4)  a  characterization  of  the  existing  pressure  propagation 
distributed  control  methodology;  (5)  suggestions  for  incorporation  into  the  existing 
control  methodology;  (6)  a  summary  of  potential  concepts  and  approaches  for  hardware, 
software  and  distributed  control  methodology  development. 

3.4  ANALYTICAL  MODEL  DEVELOPMENT 

The  models  considered  include  the  non-dynamic  and  lumped  dynamic  cases.  The  efficacy 
of  extending  the  models  to  the  more  complex  spatially  distributed  parameter  dynamic 
forms  is  contingent  on  experimental  validation  results  for  the  non-dynamic  and  lumped 
parameter  dynamic  models. 

3.4.1  Analytical  Non-Dynamic  Propagation  Delay  Model 

A  fundamental  wave  propagation  model  is  developed  for  the  case  of  catastrophic  pipe 
leakage  failure.  Consider  an  UAV  airline  of  diameter  d,  meters  that  is  fed  from  a 

compressor  incorporating  a  series  combination  of  a  debris  filter  and  a  dehydrator  baffle. 
The  effective  restrictive  diameter  of  the  filterfbaffle  combination  is  denoted  dfb  and  it  is 

assumed  that  d/b  «  dlp .  Furthermore,  let  the  UAV  airline  incur  a  large  cross-sectional 
damage  failure  of  area  af  »  nd,p  square  meters.  Under  the  prescribed  failure  conditions, 
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a  simplified  plug-flow  pressure  wavefront  is  assumed  to  characterize  the  depressurization 
propagation  process  through  the  pipeline.  The  plug-flow  process  neglects  viscous  friction 
and  assumes  that  the  delay  in  propagation,  i.e.  propagation  velocity,  is  the  result  of 
inertia.  Under  plug- flow  conditions  the  pressure  distribution  p(x,t )  at  a  point  x  »  af 
meters  away  from  the  centroid  of  the  damage  fault  may  be  expressed  the  relation 


(16) 


where  vp  is  the  plug  flow  propagation  velocity  in  meters  per  second.  The  function 
p{x,t )  may  be  expressed  in  terms  of  the  unit  step  function  u(t)  by 


p(x,t)  =  p0u(x-vpt) 


(17) 


where 


(18) 


The  propagation  of  a  pressure  waveform  through  a  piping  section  according  to  the 
inertial  propagation  delay  model  is  illustrated  in  Figure  19. 
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Figure  19:  The  propagation  of  unit  amplitude  ramped  pulse  through  a 
piping  section  with  a  dwell  time  of  six  seconds. 

3.4.2  Analytical  Lumped  Fluid  Flow  Model 

The  lumped  parameter  model  developed  incorporates  fluid  viscosity  and  compressibility 
through  Poiseuille's  law  and  the  ideal  gas  law  respectively.  The  development  considers 
the  case  of  incompressible  flow,  which  is  plausibly  applicable  for  a  wide  range  of 
commonly  encountered  air  operating  velocities. 

3.4.2. 1  Mach  Number 

Compressible  and  incompressible  fluid  flow  velocity  ranges 

The  speed  of  sound  in  a  gas  medium  is  given  by 

c  =  -JkRT  (19) 
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where  k  is  the  ratio  of  specific  heat  at  constant  pressure  to  specific  heat  at  constant 
volume  ( k-1.4  for  air),  R  is  the  universal  gas  constant  and  T  is  the  absolute  temperature. 
The  velocity  of  sound  in  air  at  STP  is  approximately  343  m/s,.  The  Mach  number  is  given 
by  the  ratio 

M  =-  (20) 

c 

where  c  is  the  speed  of  sound  and  v  is  the  flow  velocity  of  the  gas.  A  gas  flow  with  a 
Mach  number  of  0.3  or  less,  at  standard  conditions,  corresponds  to  gas  flow  velocities 
less  than  100  m/s  and  this  condition  has  been  approximated  as  incompressible  flow  [22], 
In  situations  where  the  Mach  number  is  greater  than  0.3,  the  gas  density  variation  must  be 
accounted  for. 

3. 4.2.2  Viscous  fluid  characteristics 

The  internal  friction  associated  with  a  fluid  is  characterized  by  viscosity.  This  well- 
known  property  gives  fundamental  insight  into  the  mechanism  of  the  prominent  aspects 
of  fluid  behavior  and  a  summary  of  an  approach  to  its  development  is  given,  [7],  [22], 
[25],  In  laminar  flow,  the  motion  of  a  fluid  is  considered  to  consist  of  layers  of  fluid 
sliding  past  each  other.  The  fluid  viscosity  results  in  a  frictional  coupling  between 
adjacent  layers  of  fluid.  Consequently,  a  force  must  be  exerted  between  adjacent  layers  if 
they  move  past  each  other.  One  of  the  classic  examples  of  viscous  flow  is  that  of  two 
parallel  plates  with  a  fluid  channel  separating  them.  This  situation  is  shown  in  Figure  20 
where  Plate  1  is  moving  with  velocity  v  relative  to  Plate  2  that  is  stationary.  The  fluid 
adjacent  to  Plate  2  travels  with  velocity  v  while  the  fluid  adjacent  to  Plate  2  has  velocity 
zero.  The  fluid  velocity  profile  in  the  y-direction  is  considered  to  be  distributed  among  a 
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set  of  thin  fluid  sheets,  i.e.  lamina,  which  slide  across  each  other.  The  laminas  interact 
and  transmit  force  via  the  viscosity  of  the  fluid.  In  order  to  maintain  the  motion  of  Plate  1 
a  force  must  be  continually  applied  to  this  plate.  The  force  applied  to  Plate  1  is 
transmitted  to  plate  2  through  the  viscous  interaction  of  an  infinite  series  of  laminar  sheet 
flows.  The  viscosity  therefore  is  a  measure  of  the  fluid's  resistance  to  flow. 

If  Plate  1  has  area  A,  and  the  horizontal  force  applied  to  this  plate  is  F, 
experimental  measurements  with  Newtonian  fluids  show  that  the  force  required  to 
maintain  the  velocity  of  Plate  1  is  proportional  to  the  velocity  and  plate  area,  while  being 
inversely  proportional  to  the  plate  separation  x,  [22]: 


dv  A 
dt  ~X 


(21) 


The  constant  of  proportionality  is  known  as  absolute  viscosity.  The  terms  in  equation  (21) 
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Figure  20:  Illustration  of  relative  velocity  distribution  for  parallel  plate 
laminar  flow. 
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give  rise  to  various  definitions  of  quantities  describing  fluid  physical  properties.  The 
shear  stress  exerted  on  the  fluid  is  defined  as 


Shear  stress  =  — 
A 


(22) 


When  a  shear  stress  is  applied  to  a  solid  material,  the  result  is  a  displacement  along  the 
vertical  dimension,  with  respect  to  the  configuration  of  force  shown  in  Figure  20.  Shear 
strain  is  the  ratio  of  the  resultant  displacement  dy  to  the  transverse  dimension  X  and  is 
defined  as 


Shear  strain  - 


dy 

X 


(23) 


In  the  case  of  a  fluid,  shear  strain  increases  without  limit  as  long  as  a  shear  stress  is 
applied.  In  this  case,  the  shear  stress  is  related  to  the  rate  of  change  of  the  shear  strain. 
The  rate  of  change  of  shear  strain  is  defined  as 


Strain  rate  =  — 
X 


(24) 


Furthermore,  the  viscosity  of  a  fluid  is  defined  as  the  ratio  of  shear  stress  to  strain  rate 
and  is  denoted  by 
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V 


Strain  rate 
Shear  stress 
FjA 
v/X 


(25) 


and  equation  (21)  may  be  expressed  as 


F  =  r> A  — 
A 


(26) 


The  MKS  units  of  viscosity  are  N-s-m2 .  The  practical  measure  of  viscosity  is  the  poise 
where  1  one  poise  has  MKS  units  of  (N-s-m  2)x] O' 1 .  The  viscosity  of  air  at  an  ambient 
temperature  of  20  C  is  181x1 0‘6  poise  or  1 8 1  p-poise. 

3. 4.2.3  Poiseuille’s  Law  of  viscous  flow 

As  a  viscous  fluid  flows  within  a  tube  its  velocity  varies  across  the  cross-section  of  the 
tube.  The  velocity  distribution  is  zero  at  the  wall  and  increases  to  a  maximum  value  at  the 
center  of  the  tube.  The  wall  of  the  tube  exerts  a  drag  on  the  outermost  layer  of  the  fluid 
and  this  layer  on  the  next  layer  and  so  on  and,  if  the  flow  is  laminar,  the  maximum 
velocity  occurs  along  the  center  of  the  tube.  In  effect,  this  may  be  viewed  as  a  series  of 
concentric  cylindrical  surfaces  sliding  across  each  other. 
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3.4.2.3.1  Flow  velocity 


Consider  a  pipe  of  circular  cross-section  with  radius  R  and  length  L.  Let  pressure  pi  be 
applied  to  the  left  side  of  the  pipe  and  pressure  p2  applied  to  the  right  side.  Consider  a 
cylindrical  fluid  volume  within  the  pipe  of  radius  r<R  which  is  centered  along  the  pipe 
axis,  the  net  force  F  exerted  on  the  fluid  within  this  volume  is  the  difference  between  the 
forces  applied  at  the  ends  of  this  pipe  subsection: 


F  =  —  p^Ttr1 

=  Ol  -  p2)7Tr2 


(27) 


In  steady  state,  the  velocity  profde  is  time  invariant  and  consequently  the  force 
applied  to  the  fluid  just  balances  the  viscous  retarding  force  generated  within  the  fluid  at 
the  operating  steady-state  velocity.  In  the  case  where  the  velocity  is  uniform  over  the  pipe 
cross-section,  force  is  related  to  viscosity  by  equation  (26) 

In  actuality  the  fluid  velocity  is  not  uniform  and  varies  as  a  function  of  the  radius 
r.  In  this  case,  the  relationship  between  force  and  viscosity  for  a  cylindrical  sheet  of 
length  £,  radius  r  and  thickness  dr  follows  from  (26): 

F  =  7i2nrL—.  (28) 

dr 

where  —  is  replaced  with  —  and  A  is  replaced  by  IrrrL .  Equating  the  force  in  equation 
X  dr 

(27)  to  that  of  equation  (28)  yields  the  relationship 
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(29) 


dv=  (PrJhlL 

dr  2rjL 


which  indicates  that  the  velocity  gradient,  with  respect  to  the  pipe  radius,  changes  more 
and  more  rapidly  as  the  radius  increases.  Equation  (29)  is  separable  and  characterizes  the 
velocity  profile  as 


zlAfrdr 

Ar/L  J 


(30) 


which  gives 


v ...  UhZlll 

Ar/L 


(R2-r2) 


(31) 


The  velocity  is  bounded  by  0  <  v  <  (P\P^R  with  a  maximum  value  at  the  center  of 

the  pipe  and  a  zero  velocity  at  the  pipe  wall.  The  velocity  is  proportional  to 
(P\  ~  Pi)/ ^  the  pressure  gradient,  i.e.  the  change  in  pressure  per  unit  length  of  pipe. 

3.4.2.3.2  Total  rate  of  volumetric  flow 

Consider  a  very  thin  cylindrical  surface  section  of  fluid  whose  major  axis  is  concentric 
with  that  of  a  circular  piping  line.  Let  the  inner  radius  of  the  fluid  cylinder  be  denoted  by 
r,  and  its  thickness  by  dr.  The  volume  dVvot  crossing  the  end  area  of  the  fluid  cylinder 
over  a  time  interval  dt  is  dVyoL  =  vxdAxdt  where  the  area  ,4  is  equal  to  2m-dr. 
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Substituting  equation  (31)  for  the  velocity  v  gives 


dVv0L  =—  —(R1  ~r2)27rrdrdt 


Ar/L 


(32) 


Integrating  across  the  entire  pipe  cross-section  gives  and  expression  for  the  net 
volumetric  fluid  flow  which  is  given  by 


(33) 


or 


dvvoL  _  71  \(  R*  Y  Pi  Pi  ^ 


(34) 


dt  8  rj  )y  L 


The  volumetric  flow  relationship  in  (34)  is  referred  to  as  Poiseuille's  law  [22], 

3. 4. 2. 4  Vessel  ideal  gas  equation  of  state 

Assuming  an  instantaneous  gas  fluid  equilibrium  not  subject  to  influence  of  motion, 
shear,  body  force,  stress  and  capillary  action;  the  equation  of  state  for  a  gas  is  given  by 
the  well-known  expression 


pV  =  nRT 


(35) 
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where  R  is  the  universal  gas  constant  having  a  value  of  8.314  J-motl -K\  T  is  the 
absolute  temperature,  VyoL  is  the  volume,/?  is  the  pressure,  and  n  is  the  number  of  moles 
of  the  gas  [1].  The  number  of  moles  of  gas,  n,  is  related  to  the  total  mass,  m,  of  the  gas  by 
n  =  ”/M  where  M  is  the  mass  per  mole.  Thus  (35)  may  be  expressed  in  terms  of  mass  as: 

PVvol=~RT  (36) 

M 


3.4.3  Single-ended  state  of  flow  model 

Consider  the  case  of  a  gas  vessel  being  fed  through  an  orifice  effectively  acting  as  a  feed 


Figure  21:  Pipe  section  lumped  vessel  model  subject  to  only  one  end- 
valve  open. 


pipe.  This  situation  corresponds  to  a  pipe  section  with  one  end-valve  open  and  the  other 
end-valve  closed  as  depicted  in  the  Figure  21.  The  pipe  section  is  driven  by  an  applied 
pump  pressure  p}  and  this  causes  a  fluid  mass  flow  into  the  vessel  whose  internal  fluid 
pressure  is  given  by  p2^pi-  Under  conditions  of  steady  flow,  the  continuity  of  mass  flow 
satisfies 


.  Av 

Av 

Av 

m  = - 

—  - 

— - 

V 

V 

V 

sp 

A  sp 

B  sp 

At  any  cross  section 

(37) 


50 


where  v  denotes  flow  velocity,  A  is  the  pipe  cross-sectional  area  perpendicular  to  the 
direction  of  fluid  velocity,  and  vsp  is  fluid  specific  volume.  Substituting  equation  (37)  into 
Poiseuille's  law  (34),  and  recalling  that  dVvoL  ~  vx-d.Ax.dt  yields  the  following  expression 
for  mass  flow  rate: 


dV,. 


m  = 


VOL 


dt 
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For  gases  at  a  constant  temperature 


PVsp 


A 


=  PVspB=PVsp 


Any  cross  section 


m 
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RT 


(39) 


and  furthermore,  under  constant  temperature: 


p2(  0  = 


MRT 


V, 


VES 


m(t) 


(40) 


Equations  (39)  and  (40)  may  be  respectively  expressed  as 


m(t)  =  Req (/?,  (/)  -  Pi  (t)) 


(41) 
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and 


Pi(0  =  C EQm(0 


(42) 


where  Req  represents  and  equivalent  pipe  flow  resistance  given  by 
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V  V  ) 

UJ 

(43) 


and  Ceq  represents  and  equivalent  vessel  capacitance  given  by 


C  = 
bq 


MRT 


V„ 


VES 


(44) 


3.4.3. 1  Analytical  lumped  parameter  lossless  single  end-valve  model 
Substituting  equation  (41)  into  (42)  yields  a  lumped  parameter  differential  relationship 
which  describes  a  steady-flow  pressurization  and  de-pressurization  of  the  pipe  chamber 
vessel  model: 


m{,t)  =  REQMt)-CEQm(t))  (45) 


Rearranging  (45)  gives 


m{t)  +  CEgREQm(t)  —  REqP\  ( t ) 


(46) 
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The  pressure  function  p/(t )  acts  as  an  input  or  forcing  function  to  the  pipe-vessel  system 
while  m(t )  is  the  output  or  system  response.  The  transfer  function  of  the  system  is  arrived 
at  by  taking  the  Laplace  transform  of  the  differential  equation  (46)  subject  to  zero  initial 
conditions  [2],  [4],  This  yields 


M(s) 


(47) 


The  electrical  analog  circuit  is  shown  in  Figure  22  depicts  the  charging  of  a  vessel 
due  to  a  pressure  pi(t )  applied  through  one  end  of  the  pipe. 


m(l ) 


Figure  22:  Electrical  analog  to  a  lossless  pipe  with  only  one  end- 
valve  open. 


3. 4.3. 2  Analytical  lumped  parameter  lossy  dual  end-valve  model 

The  presence  of  a  leak  precipitates  a  discharge  flow  that  reduces  the  pressure  level  within 
a  section  of  piping.  The  impact  of  a  leak  is  introduced  via  an  analogous  resistive  shunting 
section  across  the  equivalent  vessel  capacitance  as  shown  in  Figure  23.  The  piping 
section  pressure,  p2(t),  is  a  function  of  the  pressure,  pi(t),  generated  by  the  compressor 
pump(s)  and  the  leakage  rate. 
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The  piping  vessel  pressure  as  a  function  of  the  pumping  pressure  and  leakage  rate  is 
given  by  the  transfer  function 


m(0  Req 


Figure  23:  Electrical  analog  to  a  lossy  pipe  with  only  one  end 
valve  open. 
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(48) 


3.4.4  Dual-ended  state  of  flow  model 

A  transmission  pipe  within  the  system  will  generally  receive  a  supply  of  fluid  from  either 
end  as  shown  in  Figure  24.  This  may  be  the  consequence  of  the  system  being  supplied  by 
multiple  compressors  and/or  as  a  consequence  of  internal  piping  loops  present  within  the 
system.  Figure  25  illustrates  the  electrical  analog  of  the  piping  system  with  both  end- 
valves  open  in  conjunction  with  fluid  loss  piping  damage. 
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Figure  24:  Pipe  section  lumped  vessel  model  subject  to  both  end-valves  open. 

The  mass  flow  rates  and  the  vessel  pressure  are  related  to  the  compressor  driving 
pressures  and  piping  loss  by  the  Laplace  domain  matrix  equation: 


Z EQ  M  A{s) 

^ EQ  +  ^EQ _  _M B(s) _ 


Figure  25:  Electrical  analog  to  a  lossy  pipe  with  both  end-valves  open. 
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Solving  (49)  for  the  mass  flow  rates  gives 
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The  vessel  pressure  is  given  by 
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For  the  case  where  Ra=Rb=Req, 
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Solving  equation  (53)  for  the  mass  flow  rates  gives 
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Combining  equations  (54)  and  (55)  yields  and  expression  for  the  bi-directional  lossy  flow 
transfer  function 


P2(s)  =  ZEQx(MA(s)  +  MB(s)) 
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which  may  be  expressed  in  terms  of  a  multi-input  transfer  function  matrix,  H (s)f  as 


[p2(j)]  =  H(j) 
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where 
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3.4.5  Isolation  state  model  for  both  end-valves  closed 


The  remaining  response  scenario  to  be  addressed  is  the  situation  in  which  both  end  valves 
of  a  pipe  section  are  closed.  In  this  case  the  pressure  remains  constant  or,  if  a  leak  is 
present,  decays  exponentially.  The  vessel  pressure  p2(t)  is  characterized  by  the  circuit 
analog  shown  in  Figure  26. 


R 


P2(t)\  t~  initial  time  of  isolation 


-%  - 


Figure  26:  Electrical  analog  to  a  lossy  pipe  with  both  end-valves 
closed. 


From  the  instant  that  both  valves  are  closed,  the  response  of  the  piping  section  is 
independent  of  the  rest  of  the  system.  The  response  of  the  pipe  section  is  driven  by  the 
value  of  the  pipe  section  pressure  at  the  instant  that  both  valves  attained  the  closed 
position.  This  pressure  is  the  initial  condition  driving  the  homogeneous  response.  The 
homogeneous  response  is  characterized  by  the  relationship 


(59) 


58 


or  equivalently 


0  —  Pj  (0 Ceq  "*■  n  Pi  (0 
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EQ 


The  solution  to  equation  (60)  gives  the  isolated  homogeneous  pressure  response  of  the 
pipe  section: 


p2(t)  =  p2(0)xe-,/(«^) 


(61) 


3.5  SIMULATION  MODEL 

The  simulation  models  discussed  are  the  non-dynamic  propagation  delay  type  and  the 
lumped  parameter  type.  These  simulation  models  follow  from  the  analytical  models 
developed  for  the  non-dynamic  case,  equations  (16),  (17)  and  (18)  and  for  the  lumped 
dynamic  case,  equations  (48),  (58)  and  (60).  Contingent  upon  experimental  model 
validation  results,  it  may  be  necessary  to  extend  these  models  to  include  spatially 
distributed  parameters. 

3.5.1  Non-Dynamic  Propagation  Delay  Simulation  Model 

The  simulation  model  for  the  non-dynamic  case  is  shown  in  Figure  27.  The  model  fluid 
supply  is  located  on  the  left-hand  side;  the  pressure  from  this  fluid  source  is  transmitted 
to  the  right  through  two  piping  sections.  The  piping  sections  may  incorporate  the 
propagation  delay  via  a  time  delay  element.  End-valves  and  catastrophic  pipe  faults  are 
modeled  using  a  product  element  which  essentially  blocks  the  propagation  of  pressure  to 
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the  right.  This  simulation  approach  is  characterized  by  unidirectional  pressure 
propagation.  Since  the  pressure  propagation  model  employs  time  delays  a  direct 
interconnection  of  left  and  right  unidirectional  propagation  models  is  feasible  employing 
AND  gate  pressure  logic  combinatory,  since  the  time  delays  would  combine  nonlinearly. 
Micro-controller  action  to  regulate  the  end-valves  is  also  incorporated  into  the  simulation 
model.  These  consist  of  data  storage  and  logical  operator  elements  that  sample  the 
pressure  in  each  piping  section.  The  control  logic  configuration  show  will  only  allow  a 
valve  to  remain  open  if  pressure  can  be  maintained  by  both  of  the  pipe  sections  connected 
to  the  value.  While  suitable  for  the  evaluation  of  some  aspects  of  system  performance,  the 
model  lacks  the  capability  to  directly  integrate  into  an  overall  system  simulation  model 
which  incorporates  physical  fluid  phenomena  directly.  Due  to  this  fact,  an  extension  to 
dynamic  modeling  appears  warranted. 


3.5.2  Lumped  Dynamic  Simulation  Model 

Analytical  lumped  parameter  models  have  been  developed  to  describe  the  piping  vessel 
pressure  response  to  three  operational  cases.  The  first  case  corresponds  to  one  end  valve 
open  and  one  end  valve  closed  and  is  derived  from  equation  (48)  as 
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The  second  case  corresponds  to  the  condition  that  exists  when  both  end  valves  are  open. 
This  follows  from  equation  (57)  and  is  given  by 
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The  third  and  final  case  is  the  isolation  response  that  results  when  both  end-valves  are 
closed.  This  follows  from  equation  (59)  as 


0  =  sP2(s)  +  — ^ —  P2  (s)  (64) 

L'-'EQ 

Equations  (62),  (63)  and  (64)  may  be  rearranged  respectively  into  a  common  format 
suitable  for  generating  a  simulation  diagram  model  with  switched  gains: 
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P2(s)  =  P2(s) 


(67) 


where 


(68) 


and 


(69) 


The  rEQ  is  the  time  constant  associated  with  non-isolation,  i.e.  at  least  one  of  the  piping 
section  end-valves  is  open,  while  recharacterizes  the  isolation  case,  i.e.  when  both  of 
the  end  valves  are  closed. 

Figure  30  shows  the  lumped  parameter  dynamic  simulation  model.  This  model 
represents  a  single  piping  section  with  end-valves.  The  simulation  structure  integrates  the 
analytical  models  developed  for  the  pipe  section  subject  to  only  one  end-valve  open,  both 
end-valves  open  and  no  end-valves  open,  i.e.,  isolation.  The  model  is  based  on  equations 
(66),  (67),  (68).  These  equations  possess  an  identical  simulation  structure  and  the 
simulation  model  switches  between  them  based  on  the  system  valve  configurations:  no, 
only  one,  or  both  end-valves  open.  The  logic  elements  included  decide  which  model  is 
applicable  based  on  end-valve  control  signals.  In  addition  groups  of  this  pipe/valve  model 
segment  are  readily  integrated  into  a  complete  UAV  system  model.  This  entails  simply 
interconnecting  the  model  segment  pressure  inputs  to  the  pressure  output  of  adjacent  pipe 
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sections.  Furthermore,  any  number  of  piping  sections  may  be  placed  into  a  series-parallel 
topology. 

An  LP  pipeline  simulation  model  was  constructed  from  the  lumped  dynamic 
sectional  model  previously  presented.  The  pressure  response  behavior,  which 
characterizes  an  LP  air  system,  is  illustrated  in  Figures  28  and  29.  Figure  28  illustrates 
the  charging  of  a  40  foot  long  pipe  section  to  100  psi  with  one  end  of  the  pipeline  closed 
and  the  other  end  driven  by  a  compressor.  The  pressure  response  along  the  pipeline  is 
shown  for  a  two  second  interval.  In  Figure  29,  the  LP  air  system  begins  charging  as  it  did 
in  Figure  28.  At  a  time  of  one  second  the  pipeline  undergoes  a  rupture  failure  at  a 
distance  of  twenty  feet  along  the  pipe  with  a  cross  sectional  area  large  compared  to  the 
pipeline  diameter.  The  pipeline  is  observed  to  discharge  to  ambient  pressure  on  the  closed 
end  side  and  to  go  into  a  redistributed  steady  state  flow  on  the  supply  side. 


Figure  28:  Simulation  of  distributed  pressure  profile  along  a  charging  LP  air  pipeline. 
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Figure  29:  Simulation  illustrating  a  charging  LP  pipeline  that  undergoes  a 
catastrophic  rupture  failure  located  at  distance  20ft.  and  occurring  at  time  1  sec. 

3.5.3  Lumped  Dynamic  Leak  Detection  and  Characterization 

The  process  of  leak  detection  and  characterization  must  be  accomplished  with  a  set  of 
pressure  switch  sensors.  This  sensor  complement  does  not  provide  graduated  pressure 
readings  as  may  be  employed  in  a  response  fitting  algorithm  such  as  least  mean  square  or 
recursive  least  squares.  Instead,  the  leak  estimation  problem  must  be  solved  in  tenns  of 
the  switching  times  of  the  pressure  level  switch  sensors.  These  may  be  projected  to  yield 
an  effective  leakage  decay  time  constant.  The  leakage  decay  time  constant  may  be 
contrasted  with  the  baseline  non-leakage  effective  time  constant.  The  magnitude  of  the 
deviation  in  the  time  constants  yields  an  estimate  of  the  severity  of  the  leakage  failure. 
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Figure  30:  Lumped  parameter  dynamic  simulation  model. 


3.6  Development  Methodology  for  Autonomous  Damage  Mitigation 

3.6.1  Implementation  Approach 

The  fundamental  problem  of  autonomous  control  development  for  distributed  shipboard 
UAV  systems  is  next  addressed.  An  environmentally  friendly  and  cost  effective  land- 
based  UAV  system  was  chosen  as  the  initial  development  platform  with  results  extending 
directly  to  other  distributed  shipboard  systems,  including  potential  hull  structure  damage 
assessment  and  control  applications.  A  real-time  simulation/emulation  based 
methodology  is  introduced  and  employed  to  implement  a  prototype  acoustic  sensor  array 
with  the  capability  of  detecting  and  locating  acoustic  sources  within  a  distributed 
auxiliary  supply  system  and  has  the  potential  to  extend  to  more  general  and  diverse 
distributed  shipboard  systems.  The  methodology  is  based  on  simulation/emulation  and 
employs  physically  based  modeling  to  characterize  signals  and  the  signal  environment. 
The  impetus  for  developing  such  a  system  is  to  devise  a  viable  means  of  economically 
supplementing  existing  widely  used  pressure  switch  type  sensors  which  are  characterized 
by  prohibitively  long  inertial  lag  responses  precluding  their  use  for  dynamic  fault 
detection  and  location.  The  use  of  non-supplemented  switch  type  sensors  requires  quasi¬ 
steady  switching  to  isolate  faults  and  suffers  from  associated  additional  working  fluid 
loss  and  the  potential  for  extended  algorithm  convergence  time.  The  prototype 
implementation  employs  simplified  filtered  threshold  detection  for  event  recognition  and 
feature  delay  estimation.  The  filtering  provides  band  limiting  and  along  with 
environmental  interference  and  noise  reduction.  The  combination  of  reduced  bandwidth 
coupled  with  simplified  detection  is  aimed  at  reducing  the  computational  and  component 
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costs  associated  with  higher  bandwidth  detection.  Feature  delay  estimates  are  utilized  in 
sensor-to-sensor  feature  propagation  models  to  estimate  the  location  coordinates  of  a 
fault  source.  Additionally,  some  aspects  of  implementing  sensor  nodes  on  other 
distributed  shipboard  systems  and  structures  are  considered. 

3.6.2  Implementation  Background 

Aspects  of  dynamic  modeling  and  switch  sensor  based  autonomous  damage  control  for 
UAV  systems  have  been  presented  in  [29],  [10].  The  general  distributed  architecture 
considered  in  these  efforts  has  been  adopted  for  the  present  application.  A  brief 
descriptive  summary  of  the  distributed  architecture  based  on  the  above  works  provides  a 
necessary  perspective  for  the  application. 

The  layout  of  a  typical  shipboard  distributed  fluid  supply  system  configuration 
was  shown  in  Figure  18.  The  system  is  comprised  of  an  interconnection  of  series-parallel 
fluid  flow  piping  paths  with  associated  valves,  sensors,  micro-controllers  and  fluid 
supplies.  The  fluid  supply  units  are  compressors  which  feed  the  distributed  system  and 
are  assumed  to  have  relatively  high  volume  baffle  reservoirs.  The  flow  control  valve 
action  is  binary  (e.g.  on/off)  and  actuation  is  controlled  by  an  associated  micro-controller 
unit.  The  controllers  each  execute  an  independent  algorithm  whose  action  is  based  on  two 
pressure  switch  sensors  located  on  each  side  of  the  valve.  Any  of  the  piping  sections  and 
associated  valves  may  experience  a  leakage  failure.  In  this  scenario,  the  operational  goal 
is  to  autonomously  mitigate  damage  impact  by  minimizing  fluid  loss  and  maximizing 
fluid  supply  to  meet  load  demand. 
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A  simulation  model  suitable  for  the  analysis,  design  and  implementation  of 
autonomous  control  for  switch  sensor  distributed  fluid  supply  systems  has  been 
developed  in  [10].  The  simulation  is  based  on  differential  lumped  modeling  and 
incorporates  the  valve  control  actions  of  distributed  micro  controllers  and 
leakage/blockage  failure.  The  model  predicts  dynamic  distributed  UAV  pipeline  pressure 
and  is  based  on  simplified  stagnation  and  flow  equilibrium  relations.  The  model  is  readily 
modified  to  incorporate  other  types  of  fluids  and  nonlinear  flow  pressure  relationships.  In 
addition,  the  model  provides  a  basis  for  failure  identification,  condition  assessment  and 
autonomous  damage  control. 

In  the  above  UAV  context,  the  present  problem  addresses  the  dynamic  limitations 
of  pressure  switch  sensors  precluding  response  to  system  transitional  dynamics  which  are 
precipitated  by  a  sudden  catastrophic  damage  failure  in  the  piping  and/or  valves.  This 
inability  to  sense  rapid  dynamics  results  in  the  loss  information  useful  for  the  rapid 
location  and  assessment  of  system  damage.  However,  the  process  of  leak  detection 
damage  control  can  be  accomplished  with  a  set  of  pressure  switch  sensors  as  done  by 
[30]  using  quasi-steady  state  sensing.  A  more  detailed  characterization  of  the  nature  of 
pressure  switch  sensor  response  follows. 

3.6.2. 1  Exiting  Switch  Sensor  Characteristics 

The  pressure  sensors  associated  with  each  piping  valve  detect  whether  the  air  pressure 
within  the  local  section  of  pipe  is  above  or  below  a  reference  threshold  pressure  range. 
The  pressure  sensor  output  is  binary  and  provided  by  the  opening  and  closing  of  a  set  of 
switch  contacts.  The  pressure  sensor  input-output  relationship  between  local  pipe  section 
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pressure  and  the  sensor  switch  setting  exhibits  a  hysteresis  characteristic  as  illustrated  in 
Figure  3 1 .  The  sensor  threshold  pressure  is  defined  to  be 

T  +  T 

Tsr  s  a  (70) 


Figure  31:  Pressure  sensor  switch  response  as  a  function  of  quasi¬ 
steady  state  local  pipe  pressure. 


The  conservation  of  mass  is  applicable  to  both  dynamic  and  steady  state  responses  and, 
for  fluid  velocity  v  and  fluid  density/?  ,  may  be  expressed  as 


dp  _  d(vp) 
dt  dt 


(71) 


and  in  the  case  of  incompressible  flow  the  fluid  density  is  constant  yielding 
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(72) 


0  _  d£  _  S(vp) 
dt  dt 

with  the  consequence  that  velocity  of  the  fluid  is  constant  with  respect  to  position.  This 
requires  that  the  system  remain  in  steady  state  thus  not  exhibiting  any  unsteady  behavior. 
The  application  of  time  varying  forcing  flow  and  pressure  generate  a  quasi-steady  state 
response.  In  the  case  of  a  compressible  fluid  the  inertial  lag  in  pressure  switch  sensor 
dynamics  results  in  a  complete  lack  of  response  to  and  sensing  of  the  time  varying  wave 
dynamics  of  (71)  and  quasi-steady  state  forcing  dynamics  of  (72). 

3.6.23  Fluid  Flow  Characterization 

The  relations  for  motion,  continuity  and  momentum  are  employed  to  characterize  the 
compressible  fluid  flow  within  the  UAV  system.  In  the  case  of  unsteady  state  fluid  flow, 
the  unsteady  flow  occurs  when  the  fluid  motion  at  any  point  along  the  system  is  a 
function  of  time.  Initial  assumptions  about  the  fluid  flow  are  that  it  is  isentropic  and  one¬ 
dimensional.  An  isentropic  flow  is  both  adiabatic,  i.e.  one  in  which  no  heat  is  transferred 
to  or  from  the  fluid,  and  reversible  neglecting  friction  and  other  effects.  The  assumption 
of  one-dimensional  flow  neglects  variations  in  fluid  velocity  and  pressure  transverse  to 
the  primary  flow  direction. 

Given  the  inability  of  switch  sensors  to  sense  unsteady  flow,  the  location  of  leak 
damage  must  be  determined  by  steady  state  switching  of  fluid  through  the  system 
resulting  in  added  fluid  loss  and  potentially  extended  convergence  time  of  the  switching 
algorithm.  If  the  fluid  pressure  dynamics  could  be  sensed  by  switch  sensors,  the  leak 
estimation  problem  could  be  solved  more  rapidly  in  terms  of  the  dynamic  switching  times 
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of  the  pressure  level  switch  sensors.  These  may  be  projected  to  yield  an  effective  leakage 
decay  time  constant.  The  leakage  decay  time  constant  may  be  contrasted  with  the 
baseline  non-leakage  effective  time  constant.  The  magnitude  of  the  deviation  in  the  time 
constants  yields  an  estimate  of  the  severity  and  location  of  the  leakage  failure. 

The  above  transient  dynamic  approach  to  leak  detection  necessitates  the 
implementation  of  dual  wide  bandwidth  pressure  sensors  at  each  valve  location  with 
associated  increases  in  cost  and  processing  overhead.  In  contrast,  the  approach  taken  here 
is  to  supplement  the  existing  pressure  switch  sensor  with  inexpensive  low  bandwidth,  i.e. 
low  processing  overhead,  acoustic  sensors.  The  easily  mounted  supplemental  sensors  can 
be  used  to  generate  impact  damage  location  estimates  based  upon  impact  and  valve 
switching  acoustic  dynamics.  This  approach  relies  on  the  application  of  source 
localization  techniques  subject  to  the  constraints  imposed  on  vibration  propagation  in  the 
shipboard  distributed  system  under  consideration. 

3.6.3  TRANSDUCER  CONDITIONING 

The  analysis  of  the  characteristics  of  actual  sampled  acoustic  signals  reveals  the  nature  of 
the  signal  environment  and  the  requirements  of  conditioning  to  in  attempting  to  optimize 
acoustic  source  location  algorithm  performance.  In  addition  the  conditioning  function 
may  serve  in  the  capacity  of  an  anti-aliasing  filter. 

3.6.3. 1  Signal  Environment 

A  characteristic  sampled  acoustic  signal  within  the  distributed  deployment  environment 
is  shown  in  Figure  32.  The  signal  environment  was  found  to  have  a  high  degree  of  power 
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SAMPLES 

Figure  32:  Unconditioned  acoustic  sampled  acoustic  signal 


line  interference.  The  source  signal  is  quasi-periodic  with  some  nonlinear  distortion 
introduced  by  the  acoustic  transducer  generating  the  signal.  The  interference  beat  pattern 
generated  by  the  source  signal  and  power  line  components  is  also  observable.  It  was 
observed  that  the  power  line  frequency  and  associated  harmonics  exhibited  a  presence  up 
to  a  bandwidth  of  approximately  300  Hz.  Various  processing  approaches  were  evaluated 
to  mitigate  the  impact  of  power  line  frequency  interference.  These  methods  employed 
numerical  passband  filtering,  nulling,  autocorrelation  and  cross-correlation.  This 
approach  was  found  to  have  a  prohibitive  turn  around  time  due  to  the  computational 
complexity  inherent  in  these  methods.  Consequently,  a  conditioning  circuit  was 
employed. 
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The  adjusted  frequency  response  of  the  conditioning  circuit  is  shown  in  Figure  33. 


Figure  33:  Adjusted  measured  frequency  response  of  conditioning 


This  response  compensates  for  the  acoustic  element  producing  a  net  flat  passband 
characteristic.  The  lower  cutoff  frequency  must  be  chosen  to  assure  that  power  line 
primary  and  harmonic  frequencies  are  sufficiently  suppressed.  The  upper  cutoff 
frequency  serves  the  dual  function  of  reducing  interference  from  higher  frequency 
sources  and  also  providing  anti-aliasing  filtering  for  the  sampling  process  while  assuring 
sufficient  spatial  sampling  to  resolve  the  source  position  information. 
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3.6.3.2  PROTOTYPE 


The  prototype  employs  threshold  detection  to  determine  the  presence  of  sources  and  to 
estimate  feature  propagation  delay  amongst  the  acoustic  sensor  configuration.  The  fusion 
of  the  sensor-to-sensor  feature  delay  lends  itself  to  a  graphical  representation  which  could 
be  utilized  to  generate  a  graphic  display  showing  the  deployment  environment  and 
estimated  source  location. 

3. 6.3.2. 1  Detection  and  estimation 

Characteristic  sample  signals  of  an  acoustic  burst  within  the  deployment  environment  are 
shown  in  Figure  34.  The  signals  constitute  a  signal  power  and  noise  power  within  the 
bandwidth  of  the  system.  The  objective  is  to  perform  source  detection  and  localization 
capable  of  operating  at  minimal  signal  to  noise  ratio  (SNR).  Assuming  a  Gaussian  white 
noise  process,  the  optimum  processing  approach  may  be  extrapolated  for  correlation  and 
square-law  based  detection  and  estimation  as  [9]  and  [20], 

The  propagation  delay  of  acoustic  signals  is  core  to  estimating  the  location  of 
sources.  The  propagation  velocity  of  sound  is  a  function  of  temperature  and  material,  e.g. 
for  air  the  propagation  velocity  may  be  expressed  as 
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If  the  signal  sources  are  from  a  set  of  characterized  signal  sources,  the  optimum 
processing  system  may  be  based  on  cross  correlation.  The  correlation  process  compares  a 
reference  signal,  r(t),  with  the  sampled  signal,  s(t),  employing  the  cross  correlation 
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Figure  34:  Characteristic  acoustic  burst  sample 


The  cross  correlation  will  be  a  maximum  when  T  represents  the  physical  time 
delay  due  to  propagation  of  the  acoustic  wave.  The  computational  effort  associated  with 
the  cross  correlation  function  may  be  prohibitive  in  achieving  real-time  sampling  and 
detection  turn  around  time  in  implementation. 

If  the  characteristics  of  the  acoustic  source  are  unknown  or  the  computation 
associated  with  the  cross  correlation  function  is  of  prohibitive  complexity,  square  law 
based  detection  and  estimation  may  be  a  viable  alternative.  The  processing  of  a  signal 
occurs  over  a  time  interval  r  and  over  a  bandwidth  w.  The  output  of  the  processing  may 
be  the  time  average  of  the  detected  signal  voltage  or  power  (i.e.  square  law  based 
processing).  For  the  interval  t,  a  decision  must  be  made  as  to  whether  an  acoustic  signal 
is  present.  This  is  accomplished  by  setting  a  threshold  level  for  the  output  of  the 
processing  function  over  the  processing  time  interval.  It  is  desired  to  have  the  processing 
time  interval  as  short  as  possible  for  real-time  implementation  feasibility.  The  detection 
and  performance  thresholds  for  the  cross  correlation  and  square  law  threshold  based 
approaches  are  given  by  [9]  and  [20], 
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D cross  ccrda, ion  =101og{,SAffi} 


(75) 


and 

0,,„,,,„,=51og{SVS}  (76) 

Square  law  based  detection  and  estimation  was  implemented  in  the  prototype 
since  the  acoustic  sources  are  not  characterized  sufficiently  for  cross  correlation 
computation  and  also  to  avoid  the  computational  complexity  of  cross  correlation  in  real¬ 
time  implementation.  In  the  prototype  case,  the  square  law  process  is  approximated  by 
absolute  magnitude  possessing  with  respect  to  a  threshold.  The  fundamental  problem  is  to 
trade  off  as  little  performance  as  possible  in  transitioning  from  correlation-based 
detection  and  estimation  to  that  of  square  law  based  processes. 

3.63.2.2  DOA/TDOA  Localization 

The  hyper- geometric  characterization  of  source  location  with  respect  to  a  constellation  of 
sensors  is  illustrated  in  Figure  35  without  loss  of  generality,  the  source  is  located  at 
coordinates  (8,10)  while  the  sensors  are  symmetrically  located  at  (0,-5)  and  (0,5) 
respectively  without  loss  of  generality.  The  source  is  in  a  near  field  configuration, 
approaching  the  asymptotic  region,  with  respect  to  the  sensors  [9],  [4].  The  resulting 
DOA  estimate  for  this  configuration  is  depicted  as  a  straight  line  plot  emanating  from  the 
origin.  The  intersection  of  this  plot  with  the  hyper-geometric  curve  represents  an  estimate 
of  the  source  location.  The  appropriate  hyperbolic  curve  is  selected  based  on  the  time  of 
arrival  of  the  acoustic  wave  at  the  various  sensors.  In  the  case  illustrated,  the  signal 
arrives  at  Sensor  1  first  indicating  that  the  right-side  of  the  hyperbolic  curve  contains  the 
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source.  Due  to  the  geometric  symmetry,  an  alias  source  is  present  at  coordinates  (8,-10). 
The  correct  source  location  must  be  determined  by  corroborating  TDOA  location  data 
from  other  sensor  platforms. 

The  hyper-geometric  curves  are  characterized  by  the  relation 


where  ±  c  are  the  sensor  locations,  i.e.  foci,  and  ±a  are  the  vertices  of  the  hyperbola  and 
represent  the  intersection  points  of  the  hyperbola  with  the  transverse  axis.  If  tn  is  the 
waveform  time  delay  from  sensor  1  to  sensor  2,  then 


where  v  is  the  velocity  of  wave  propagation.  If  three  or  more  sensors  are  present  then, 
given  the  relative  time  delays  associated  with  each  sensor  pair,  the  source  location  may 
be  uniquely  selected  from  the  intersection  of  two  hyperbolic  curves  given  by 


(79) 


(80) 
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The  above  approach  to  source  location  employs  a  blend  of  DOA  and  TDOA 
methods.  The  method  produces  conjugate  aliased  source  location  estimates.  Two  source 
regions  have  been  provisionally  recognized:  (1)  a  near  field  region  exhibiting  reduced 
sensitivity  and  higher  resolution;  (2)  an  asymptotic  region  (far-field)  exhibiting  greater 
sensitivity  and  reduced  resolution.  The  source  location  data  produced  is  appropriate  for 
transmission  in  a  UAV  communications  network  for  corroboration  to  remove  conjugate 
alias  source  locations. 


Figure  35:  Illustration  of  the  combined  localization  approach. 
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4.0  APPENDIX  A:  Matlab  Shell  Code 


Matlab/Mfile  {or  navigation  around  a  circular  obstacle 


clear 

simin=[0,0] 

siminl  =[0,0] 

simin2=[0,2] 

simin3=[0,2] 

path1=simin2(2) 

path2=simin3(2) 

t(1)=0 

for  i=1:1:40 

sim('circ3’) 

[r,c]=size(simout) 

[r1,c1]=size(simout1) 

[r3,c3]=size(simout3) 

[r2,c2]=size(simout2) 

[rt,ct]=size(tout) 

pathl  (i+1  )=simout2(r2) 

path2(i+1  )=simout3(r3) 

t(i+1)=t(i)+tout(rt) 

simin=[0,simout(r)] 

siminl  =[0,simout1(r1)] 

simin2=[0,simout2(r2)] 

simin3=[0,simout3(r3)j 

plot(path1,path2,’-*’) 
center=[4  4] 
radius=[.58] 

NOP=[1000] 

THETA=linspace(0,2*pi,NOP); 
RHO=ones(1  ,NOP)*radius; 

[X,  Y]  =  pol2cart(THETA,RHO); 
X=X+center(1); 

Y=Y+center(2); 

H=plot(X,Yi’r'); 
axis  square; 
hold  on 

axis([-1 ,11,-1,11]) 


grid 

end 
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4.1  APPENDIX  B  :  Overview  of  Fuzzy  Fusion  Methodology 
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